これは,円周率を巨大な桁数で計算するパッケージです. これを作ったきっかけは,ある研究で作った FFT ベースの多倍長計算 ルーチンのベンチマークを行ったことです. 計算速度は Super_PI ver 1.1 @東大金田研究室 と比較して約 2.5 倍高速です. 多倍長基本演算ルーチンは四則演算と平方根です. プログラムの変更で円周率以外の計算(sqrt(2) など)もできると思います.
Makefile : make file - for 32bit CPU: digits <= 6,000,000,000 Makefile_64bit : make file - for 64bit CPU Makefile_quad : make file - for 128bit FPU config.h : machine depending macros fftsgx.c : FFT Package in C - split-radix - use work areas fftsg_hx.c : FFT Package in C - split-radix - use no work areas pi_fftca.c : PI Calc. Program - standard version - use fft*gx.c" pi_fftcs.c : PI Calc. Program - mem save version - use "fft*g_hx.c" pi_fftcw.c : PI Calc. Program - mem swap version - use "fft*g_hx.c" dgt_div.c : data converter - from pi.dat to Super_PI format README.TXT : readme file win32bin/ : DIR win32bin/Makefile : make file - for intel C++ compiler win32bin/pi_cs.exe : Win32 exec. - mem save version win32bin/pi_cs_thread.exe : Win32 exec. - use FFT level threads win32bin/pi_cw.exe : Win32 exec. - mem swap version win32bin/pi_cw_thread.exe : Win32 exec. - use FFT level threads win32bin/dgt_div.exe : data converter - pi.dat => pi.txt a printout format
ファイル "config.h" 中にあるマクロをチェックし 必要ならば修正してください.とくに,FFT_ERROR_MARGIN は 重要なパラメータです.もし FFT_ERROR_MARGIN が非常に小さいならば 計算効率は悪くなります.もし FFT_ERROR_MARGIN >= 0.5 と設定すると 計算を間違えるかもしれません.具体的なトラブルは多倍長乗算を間違えて ニュートン反復に失敗します.FFT_ERROR_MARGIN が 10 程度でも 運がよければ計算できます.この場合, 必ず検証のための検算を行ってください.
1,048,576 桁 |
16,777,216 桁 |
134,217,728 桁 |
|
---|---|---|---|
pi_cas5 (pi_fftca.c + fftsg.c) |
54 秒 | 1323 秒 (22 分) | - |
pi_css5 (pi_fftcs.c + fftsg_h.c) |
48 秒 | 1189 秒 (19 分) | - |
pi_cws5 (pi_fftcw.c + fftsg_h.c) |
71 秒 total [65 秒 cpu] |
1702 秒 (28 分) total [1563 秒 (26 分) cpu] |
23533 秒 (6 時間 32 分) total [15748 秒 (4 時間 22 分) cpu] |
Super_PI@Kanada Lab. | 237 秒 total | 5440 秒 (90 分) total | - |
ちなみに金田,田村氏が1987年に1億3千万桁の世界記録を樹立した時の 計算時間は,当時のNECのスーパーコンピュータ SX-2 で35時間15分です. Pentium III 550 の処理速度(倍精度FLOPS)は SX-2 より遅いにもかかわらず, このプログラムは当時の彼らの5倍以上の計算速度を出しています.
1,048,576 桁 |
16,777,216 桁 |
|
---|---|---|
pi_cas5 (pi_fftca.c + fftsg.c) |
77 秒 | 2406 秒 (40 分) |
pi_css5 (pi_fftcs.c + fftsg_h.c) |
81 秒 | 2505 秒 (41 分) |
---- 計算桁数と FFT の長さとの関係 ----
ndigit = nfft*log_10(R), R >= 10000 or 1000
R は多倍長計算の進数です. R は DBL_ERROR_MARGIN と FFT,計算機,コンパイラの精度に依存します.
よく知られているように N 桁の乗算計算は筆算の方法でプログラムすると N^2 回の基本演算が必要になり,非常に計算が遅くなります.しかし, FFT を利用して計算すると N*log(N) 回ですみます. これは具体的には 1GHz の 1 CPU で10億桁の乗算を行った場合, 前者が数十年かかるのに対して,後者は数分で済みます. FFT による a と b の乗算手順は以下のとおりです.
ただし,FFT による畳み込みは巡回畳み込みなので,巡回しないように ゼロをつめるなどの処理が必要です.
- a,b の各桁を数列 a_j,b_j とおいて,それに対する 離散フーリエ変換 A_k,B_k を FFT で計算する.
- C_k = A_k * B_k を計算して C_k に対して逆 FFT を行って 畳み込み c_j を得る.
- c_j を c_j < 計算の進数 となるように正規化して 乗算 c を得る.
FFT の乗算計算の高速化に関してはさまざまなテクニックがあります. 詳しくは以下の文献等を参照してください.
円周率を計算するアルゴリズムは,大きく分けて二通りの方法があります. ひとつは,無限級数を計算する方法で,有名なものにアークタンジェント公式が あります.もうひとつは,漸化式による反復公式で,AGM (算術幾何平均法) が 有名です.
ここで使っている AGM アルゴリズムは,ガウスの算術幾何平均による 完全楕円積分の計算アルゴリズムと,楕円積分に関するルジャンドルの関係式 を組み合わせる方法 (ガウスルジャンドル公式) をさらに改良したものです. 詳しい導出および証明は 2004年数理解析研究所公開講座資料(pdf 253KB) を見て下さい.
---- a formula based on the AGM (Arithmetic-Geometric Mean) ---- c = sqrt(0.125); a = 1 + 3 * c; b = sqrt(a); e = b - 0.625; b = 2 * b; c = e - c; a = a + e; npow = 4; do { npow = 2 * npow; e = (a + b) / 2; b = sqrt(a * b); e = e - b; b = 2 * b; c = c - e; a = e + b; } while (e > SQRT_SQRT_EPSILON); e = e * e / 4; a = a + b; pi = (a * a - e - e / 2) / (a * c - e) / npow; ---- modification ---- This is a modified version of Gauss-Legendre formula (by T.Ooura). It is faster than original version.
この公式は,二次収束し一回の反復で桁数が倍になり,わずか 26回のループの反復で14億桁以上の結果が得られます.